using System;
using System.Collections.Generic;
using System.Text;
using System.IO;

namespace CddSharp
{
    public partial class Cdd
    {
        /* cddio.c:  Basic Input and Output Procedures for cddlib
           written by Komei Fukuda, fukuda@ifor.math.ethz.ch
           Version 0.94, Aug. 4, 2005
        */

        /* cddlib : C-library of the double description method for
           computing all vertices and extreme rays of the polyhedron 
           P= {x :  b - A x >= 0}.  
           Please read COPYING (GNU General Public Licence) and
           the manual cddlibman.tex for detail.
        */

#if nemkell
void SetInputFile(FILE **f,DataFileType inputfile,ErrorType *Error)
{
  int opened=0,stop,quit=0;
  int i,dotpos=0,trial=0;
  char ch;
  char *tempname;
  
  
  *Error=NoError;
  while (!opened && !quit) {
    fprintf(stderr,"\n>> Input file: ");
    scanf("%s",inputfile);
    ch=getchar();
    stop=false;
    for (i=0; i<filenamelen && !stop; i++){
      ch=inputfile[i];
      switch (ch) {
        case '.': 
          dotpos=i+1;
          break;
        case ';':  case ' ':  case '\0':  case '\n':  case '\t':     
          stop=true;
          tempname=(char*)calloc(filenamelen,sizeof(ch));
          strncpy(tempname,inputfile,i);
          strcpy(inputfile,tempname);
          free(tempname);
          break;
      }
    }
    if ( ( *f = fopen(inputfile,"r") )!= null) {
      fprintf(stderr,"input file %s is open\n",inputfile);
      opened=1;
      *Error=NoError;
    }
    else{
      fprintf(stderr,"The file %s not found\n",inputfile);
      trial++;
      if (trial>5) {
        *Error=IFileNotFound;
        quit=1;
      }
    }
  }
}

void SetWriteFileName(DataFileType inputfile, DataFileType outfile, char cflag, RepresentationType rep)
{
  char *extension;
  DataFileType ifilehead="";
  int i,dotpos;
  
  switch (cflag) {
    case 'o':
      switch (rep) {
        case Generator:
          extension=".ine"; break;     /* output file for ine data */
        case Inequality:
          extension=".ext"; break;     /* output file for ext data */
        default: 
          extension=".xxx";break;
      }
      break;

    case 'a':         /* decide for output adjacence */
      if (rep==Inequality)
        extension=".ead";       /* adjacency file for ext data */
      else
        extension=".iad";       /* adjacency file for ine data */
      break;
    case 'i':         /* decide for output incidence */
      if (rep==Inequality)
        extension=".ecd";       /* ext incidence file */
      else
        extension=".icd";       /* ine incidence file */
      break;
    case 'n':         /* decide for input incidence */
      if (rep==Inequality)
        extension=".icd";       /* ine incidence file */
      else
        extension=".ecd";       /* ext incidence file */
      break;
    case 'j':        /* decide for input adjacency */
      if (rep==Inequality)
        extension=".iad";       /* ine adjacency file */
      else
        extension=".ead";       /* ext adjacency file */
      break;
    case 'l':
      extension=".ddl";break;   /* log file */
    case 'd':
      extension=".dex";break;   /* decomposition output */
    case 'p':
      extension="sub.ine";break;  /* preprojection sub inequality file */
    case 'v':
      extension=".solved";break;  /* verify_input file */
    case 's':
      extension=".lps";break;   /* LP solution file */
    default:
      extension=".xxx";break;
  }
  dotpos=-1;
  for (i=0; i< strlen(inputfile); i++){
    if (inputfile[i]=='.') dotpos=i;
  }
  if (dotpos>1) strncpy(ifilehead, inputfile, dotpos);
  else strcpy(ifilehead,inputfile);
  if (strlen(inputfile)<=0) strcpy(ifilehead,"tempcdd");
  strcpy(outfile,ifilehead); 
  strcat(outfile,extension); 
  if (strcmp(inputfile, outfile)==0) {
    strcpy(outfile,inputfile); 
    strcat(outfile,extension); 
  }
/*  fprintf(stderr,"outfile name = %s\n",outfile);  */
}


NumberType GetNumberType(char *line)
{
  NumberType nt;

  if (strncmp(line, "integer", 7)==0) {
    nt = Integer;
  }
  else if (strncmp(line, "rational", 8)==0) {
    nt = Rational;
  }
  else if (strncmp(line, "real", 4)==0) {
    nt = Real;
  }
  else { 
    nt=Unknown;
  }
  return nt;
}

void ProcessCommandLine(FILE *f, Matrix M, char *line)
{
  char newline[linelenmax];
  int j;
  double value;

  init(value);
  if (strncmp(line, "hull", 4)==0) {
    M.representation = Generator;
  }
  if (strncmp(line, "debug", 5)==0) {
    debug = true;
//#ifdef GMPRATIONAL
    ddf_debug = ddf_TRUE;
//#endif
  }
  if (strncmp(line, "partial_enum", 12)==0 ||
       strncmp(line, "equality", 8)==0  ||
       strncmp(line, "linearity", 9)==0 ) {
    fgets(newline,linelenmax,f);    
    SetLinearity(M,newline);
  }
  if (strncmp(line, "maximize", 8)==0 ||
      strncmp(line, "minimize", 8)==0) {
    if (strncmp(line, "maximize", 8)==0) M.objective=LPmax;
    else M.objective=LPmin;
    for (j = 1; j <= M.colsize; j++) {
    if (M.numbtype==Real) {
//#if !defined(GMPRATIONAL)
        double rvalue;
        fscanf(f, "%lf", &rvalue);
        set_d(value, rvalue);
//#endif
      } else {
        fread_rational_value (f, value);
      }
      set(M.rowvec[j - 1],value);
      if (debug) {fprintf(stderr,"cost(%5ld) =",j); WriteNumber(stderr,value);}
    }  /*of j*/
  }
  clear(value);
}
#endif

        public bool AppendMatrix2Poly(Polyhedra poly, Matrix M)
        {
            bool success = false;
            Matrix Mpoly, Mnew = null;
            ErrorType err;

            if (poly != null && poly.m >= 0 && poly.d >= 0 &&
                poly.d == M.colsize && M.rowsize > 0)
            {
                Mpoly = CopyInput(poly);
                Mnew = AppendMatrix(Mpoly, M);

                poly = DDMatrix2Poly(Mnew, out err);

                if (err == ErrorType.NoError) success = true;
            }
            return success;
        }

        public Matrix MatrixCopy(Matrix M)
        {
            Matrix Mcopy = null;
            int m;
            int d;

            m = M.rowsize;
            d = M.colsize;
            if (m >= 0 && d >= 0)
            {
                Mcopy = CreateMatrix(m, d);
                CopyAmatrix(Mcopy.matrix, M.matrix, m, d);
                CopyArow(Mcopy.rowvec, M.rowvec, d);
                M.linset.CopyTo(Mcopy.linset);
                Mcopy.numbtype = M.numbtype;
                Mcopy.representation = M.representation;
                Mcopy.objective = M.objective;
            }
            return Mcopy;
        }

        public Matrix CopyMatrix(Matrix M)
        {
            return MatrixCopy(M);
        }

        public Matrix MatrixNormalizedCopy(Matrix M)
        {
            Matrix Mcopy = null;
            int m;
            int d;

            m = M.rowsize;
            d = M.colsize;
            if (m >= 0 && d >= 0)
            {
                Mcopy = CreateMatrix(m, d);
                CopyNormalizedAmatrix(Mcopy.matrix, M.matrix, m, d);
                CopyArow(Mcopy.rowvec, M.rowvec, d);
                M.linset.CopyTo(Mcopy.linset);
                Mcopy.numbtype = M.numbtype;
                Mcopy.representation = M.representation;
                Mcopy.objective = M.objective;
            }
            return Mcopy;
        }


        public Matrix MatrixAppend(Matrix M1, Matrix M2)
        {
            Matrix M = null;
            int i, m, m1, m2;
            int j, d, d1, d2;

            m1 = M1.rowsize;
            d1 = M1.colsize;
            m2 = M2.rowsize;
            d2 = M2.colsize;

            m = m1 + m2;
            d = d1;

            if (d1 >= 0 && d1 == d2 && m1 >= 0 && m2 >= 0)
            {
                M = CreateMatrix(m, d);
                CopyAmatrix(M.matrix, M1.matrix, m1, d);
                CopyArow(M.rowvec, M1.rowvec, d);
                for (i = 0; i < m1; i++)
                {
                    if (M1.linset.Get(i + 1)) M.linset.Set(i + 1, true);
                }
                for (i = 0; i < m2; i++)
                {
                    for (j = 0; j < d; j++)
                        M.matrix[m1 + i][j] = M2.matrix[i][j];
                    /* append the second matrix */
                    if (M2.linset.Get(i + 1)) M.linset.Set(m1 + i + 1, true);
                }
                M.numbtype = M1.numbtype;
            }
            return M;
        }

        public Matrix MatrixNormalizedSortedCopy(Matrix M, out int[] newpos)  /* 094 */
        {
            /* Sort the rows of double[][] lexicographically, and return a link to this sorted copy. 
            The vector newpos is allocated, where newpos[i] returns the new row index
            of the original row i (i=1,...,M.rowsize). */
            Matrix Mcopy = null, Mnorm = null;
            int m, i;
            int d;
            int[] roworder;

            /* if (newpos!=null) free(newpos); */
            m = M.rowsize;
            d = M.colsize;
            roworder = new int[m + 1];
            newpos = new int[m + 1];
            if (m >= 0 && d >= 0)
            {
                Mnorm = MatrixNormalizedCopy(M);
                Mcopy = CreateMatrix(m, d);
                for (i = 1; i <= m; i++) roworder[i] = i;

                RandomPermutation(roworder, m, 123);
                QuickSort(roworder, 1, m, Mnorm.matrix, d);

                PermuteCopyAmatrix(Mcopy.matrix, Mnorm.matrix, m, d, roworder);
                CopyArow(Mcopy.rowvec, M.rowvec, d);
                for (i = 1; i <= m; i++)
                {
                    if (M.linset.Get(roworder[i])) Mcopy.linset.Set(i, true);
                    newpos[roworder[i]] = i;
                }
                Mcopy.numbtype = M.numbtype;
                Mcopy.representation = M.representation;
                Mcopy.objective = M.objective;

            }

            return Mcopy;
        }

        public Matrix MatrixUniqueCopy(Matrix M, out int[] newpos)
        {
            /* Remove row duplicates, and return a link to this sorted copy. 
               Linearity rows have priority over the other rows.
               It is better to call this after sorting with MatrixNormalizedSortedCopy.
               The vector newpos is allocated, where *newpos[i] returns the new row index
               of the original row i (i=1,...,M.rowsize).  *newpos[i] is negative if the original
               row is dominated by -*newpos[i] and eliminated in the new copy.
            */
            Matrix Mcopy = null;
            int m, i, uniqrows;
            BitArray preferredrows;
            int d;
            int[] roworder;

            /* if (newpos!=null) free(newpos); */
            m = M.rowsize;
            d = M.colsize;
            preferredrows = M.linset;
            roworder = new int[m + 1];
            if (m >= 0 && d >= 0)
            {
                for (i = 1; i <= m; i++) roworder[i] = i;
                UniqueRows(roworder, 1, m, M.matrix, d, preferredrows, out uniqrows);

                Mcopy = CreateMatrix(uniqrows, d);
                PermutePartialCopyAmatrix(Mcopy.matrix, M.matrix, m, d, roworder, 1, m);
                CopyArow(Mcopy.rowvec, M.rowvec, d);
                for (i = 1; i <= m; i++)
                {
                    if (roworder[i] > 0 && M.linset.Get(i)) Mcopy.linset.Set(roworder[i], true);
                }
                Mcopy.numbtype = M.numbtype;
                Mcopy.representation = M.representation;
                Mcopy.objective = M.objective;
            }
            newpos = roworder;
            return Mcopy;
        }


        public Matrix MatrixNormalizedSortedUniqueCopy(Matrix M, out int[] newpos) /* 094 */
        {
            /* Sort and remove row duplicates, and return a link to this sorted copy. 
               Linearity rows have priority over the other rows.
               It is better to call this after sorting with MatrixNormalizedSortedCopy.
               The vector newpos is allocated, where *newpos[i] returns the new row index
               of the original row i (i=1,...,M.rowsize).  *newpos[i] is negative if the original
               row is dominated by -*newpos[i] and eliminated in the new copy.
            */
            Matrix M1 = null, M2 = null;
            int m, i;
            int d;
            int[] newpos1 = null, newpos1r = null, newpos2 = null;

            /* if (newpos!=null) free(newpos); */
            m = M.rowsize;
            d = M.colsize;
            newpos = new int[m + 1];
            newpos1r = new int[m + 1];
            if (m >= 0 && d >= 0)
            {
                M1 = MatrixNormalizedSortedCopy(M, out newpos1);
                for (i = 1; i <= m; i++) newpos1r[newpos1[i]] = i;  /* reverse of newpos1 */
                M2 = MatrixUniqueCopy(M1, out newpos2);
                M2.linset.SetAll(false);    //emptyset
                for (i = 1; i <= m; i++)
                {
                    if (newpos2[newpos1[i]] > 0)
                    {
                        //printf("newpos1[%ld]=%ld, newpos2[newpos1[%ld]]=%ld\n",i,newpos1[i], i,newpos2[newpos1[i]]);
                        if (M.linset.Get(i)) M2.linset.Set(newpos2[newpos1[i]], true);
                        newpos[i] = newpos2[newpos1[i]];
                    }
                    else
                    {
                        newpos[i] = -newpos1r[-newpos2[newpos1[i]]];
                    }
                }

            }

            return M2;
        }

        public Matrix MatrixSortedUniqueCopy(Matrix M, out int[] newpos)  /* 094 */
        {
            /* Same as MatrixNormalizedSortedUniqueCopy except that it returns a unnormalized origial data
               with original ordering.
            */
            Matrix M1 = null, M2 = null;
            int m, i, k, ii;
            int d;
            int[] newpos1 = null, newpos1r = null, newpos2 = null;

            /* if (newpos!=null) free(newpos); */
            m = M.rowsize;
            d = M.colsize;
            newpos = new int[m + 1];
            newpos1r = new int[m + 1];
            if (m >= 0 && d >= 0)
            {
                M1 = MatrixNormalizedSortedCopy(M, out newpos1);
                for (i = 1; i <= m; i++) newpos1r[newpos1[i]] = i;  /* reverse of newpos1 */
                M2 = MatrixUniqueCopy(M1, out newpos2);

                M2.linset.SetAll(false);
                for (i = 1; i <= m; i++)
                {
                    if (newpos2[newpos1[i]] > 0)
                    {
                        if (M.linset.Get(i)) M2.linset.Set(newpos2[newpos1[i]], true);
                        newpos[i] = newpos2[newpos1[i]];
                    }
                    else
                    {
                        newpos[i] = -newpos1r[-newpos2[newpos1[i]]];
                    }
                }

                ii = 0;
                M2.linset.SetAll(false);
                for (i = 1; i <= m; i++)
                {
                    k = newpos[i];
                    if (k > 0)
                    {
                        ii += 1;
                        newpos[i] = ii;
                        CopyArow(M2.matrix[ii - 1], M.matrix[i - 1], d);
                        if (M.linset.Get(i)) M2.linset.Set(ii, true);
                    }
                }

            }

            return M2;
        }

        public Matrix AppendMatrix(Matrix M1, Matrix M2)
        {
            return MatrixAppend(M1, M2);
        }

        public bool MatrixAppendTo(ref Matrix M1, Matrix M2)
        {
            Matrix M = null;
            int i, m, m1, m2;
            int j, d, d1, d2;
            bool success = false;

            m1 = M1.rowsize;
            d1 = M1.colsize;
            m2 = M2.rowsize;
            d2 = M2.colsize;

            m = m1 + m2;
            d = d1;

            if (d1 >= 0 && d1 == d2 && m1 >= 0 && m2 >= 0)
            {
                M = CreateMatrix(m, d);
                CopyAmatrix(M.matrix, M1.matrix, m1, d);
                CopyArow(M.rowvec, M1.rowvec, d);
                for (i = 0; i < m1; i++)
                {
                    if (M1.linset.Get(i + 1)) M.linset.Set(i + 1, true);
                }
                for (i = 0; i < m2; i++)
                {
                    for (j = 0; j < d; j++)
                        M.matrix[m1 + i][j] = M2.matrix[i][j];
                    /* append the second matrix */
                    if (M2.linset.Get(i + 1)) M.linset.Set(m1 + i + 1, true);
                }
                M.numbtype = M1.numbtype;

                M1 = M;
                success = true;
            }
            return success;
        }

        public bool MatrixRowRemove(Matrix M, int r) /* 092 */
        {
            int i, m;
            int d;
            bool success = false;

            m = M.rowsize;
            d = M.colsize;

            if (r >= 1 && r <= m)
            {
                M.rowsize = m - 1;

                M.linset.Set(r, false);
                /* slide the row headers */
                for (i = r; i < m; i++)
                {
                    M.matrix[i - 1] = M.matrix[i];
                    if (M.linset.Get(i + 1))
                    {
                        M.linset.Set(i + 1, false);
                        M.linset.Set(i, true);
                    }
                }
                success = true;
            }
            return success;
        }

        public bool MatrixRowRemove2(Matrix M, int r, out int[] newpos) /* 094 */
        {
            int i, m;
            int d;
            bool success = false;
            int[] roworder;

            newpos = null;

            m = M.rowsize;
            d = M.colsize;

            if (r >= 1 && r <= m)
            {
                roworder = new int[m + 1];
                M.rowsize = m - 1;

                M.linset.Set(r, false);
                /* slide the row headers */
                for (i = 1; i < r; i++) roworder[i] = i;
                roworder[r] = 0; /* meaning it is removed */
                for (i = r; i < m; i++)
                {
                    M.matrix[i - 1] = M.matrix[i];
                    roworder[i + 1] = i;
                    if (M.linset.Get(i + 1))
                    {
                        M.linset.Set(i + 1, false);
                        M.linset.Set(i, true);
                    }
                }
                success = true;
            }
            return success;
        }

        public Matrix MatrixSubmatrix(Matrix M, BitArray delset) /* 092 */
        {
            Matrix Msub = null;
            int i, isub = 1, m, msub;
            int d;

            m = M.rowsize;
            d = M.colsize;
            msub = m;
            if (m >= 0 && d >= 0)
            {
                for (i = 1; i <= m; i++)
                {
                    if (delset.Get(i)) msub -= 1;
                }
                Msub = CreateMatrix(msub, d);
                for (i = 1; i <= m; i++)
                {
                    if (!delset.Get(i))
                    {
                        CopyArow(Msub.matrix[isub - 1], M.matrix[i - 1], d);
                        if (M.linset.Get(i))
                        {
                            Msub.linset.Set(isub, true);
                        }
                        isub++;
                    }
                }
                CopyArow(Msub.rowvec, M.rowvec, d);
                Msub.numbtype = M.numbtype;
                Msub.representation = M.representation;
                Msub.objective = M.objective;
            }
            return Msub;
        }

        public Matrix MatrixSubmatrix2(Matrix M, BitArray delset, out int[] newpos) /* 092 */
        { /* returns a pointer to a new matrix which is a submatrix of M with rows in delset
  removed.  *newpos[i] returns the position of the original row i in the new matrix.
  It is -1 if and only if it is deleted.
  */

            Matrix Msub = null;
            int i, isub = 1, m, msub;
            int d;
            int[] roworder;

            newpos = null;

            m = M.rowsize;
            d = M.colsize;
            msub = m;
            if (m >= 0 && d >= 0)
            {
                roworder = new int[m + 1];
                for (i = 1; i <= m; i++)
                {
                    if (delset.Get(i)) msub -= 1;
                }
                Msub = CreateMatrix(msub, d);
                for (i = 1; i <= m; i++)
                {
                    if (delset.Get(i))
                    {
                        roworder[i] = 0; /* zero means the row i is removed */
                    }
                    else
                    {
                        CopyArow(Msub.matrix[isub - 1], M.matrix[i - 1], d);
                        if (M.linset.Get(i))
                        {
                            Msub.linset.Set(isub, true);
                        }
                        roworder[i] = isub;
                        isub++;
                    }
                }
                newpos = roworder;
                CopyArow(Msub.rowvec, M.rowvec, d);
                Msub.numbtype = M.numbtype;
                Msub.representation = M.representation;
                Msub.objective = M.objective;
            }
            return Msub;
        }


        public Matrix MatrixSubmatrix2L(Matrix M, BitArray delset, out int[] newpos) /* 094 */
        { /* This is same as MatrixSubmatrix2 except that the linearity rows will be shifted up
     so that they are at the top of the matrix.
  */
            Matrix Msub = null;
            int i, iL, iI, m, msub;
            int d;
            int[] roworder;

            newpos = null;

            m = M.rowsize;
            d = M.colsize;
            msub = m;
            if (m >= 0 && d >= 0)
            {
                roworder = new int[m + 1];
                for (i = 1; i <= m; i++)
                {
                    if (delset.Get(i)) msub -= 1;
                }
                Msub = CreateMatrix(msub, d);
                iL = 1; iI = M.linset.Count + 1;  /* starting positions */
                for (i = 1; i <= m; i++)
                {
                    if (delset.Get(i))
                    {
                        roworder[i] = 0; /* zero means the row i is removed */
                    }
                    else
                    {
                        if (M.linset.Get(i))
                        {
                            CopyArow(Msub.matrix[iL - 1], M.matrix[i - 1], d);
                            Msub.linset.Set(i, false);
                            Msub.linset.Set(iL, true);
                            roworder[i] = iL;
                            iL += 1;
                        }
                        else
                        {
                            CopyArow(Msub.matrix[iI - 1], M.matrix[i - 1], d);
                            roworder[i] = iI;
                            iI += 1;
                        }
                    }
                }
                newpos = roworder;
                CopyArow(Msub.rowvec, M.rowvec, d);
                Msub.numbtype = M.numbtype;
                Msub.representation = M.representation;
                Msub.objective = M.objective;
            }
            return Msub;
        }

        public int MatrixRowsRemove(ref Matrix M, BitArray delset) /* 094 */
        {
            Matrix Msub = null;
            int success;

            Msub = MatrixSubmatrix(M, delset);

            M = Msub;
            success = 1;
            return success;
        }

        public int MatrixRowsRemove2(ref Matrix M, BitArray delset, out int[] newpos) /* 094 */
        {
            Matrix Msub = null;
            int success;

            Msub = MatrixSubmatrix2(M, delset, out newpos);

            M = Msub;
            success = 1;
            return success;
        }

        public int MatrixShiftupLinearity(ref Matrix M, out int[] newpos) /* 094 */
        {
            Matrix Msub = null;
            int success;
            BitArray delset;

            delset = new BitArray(M.rowsize);  /* emptyset */
            Msub = MatrixSubmatrix2L(M, delset, out newpos);

            M = Msub;

            success = 1;
            return success;
        }

        public Polyhedra CreatePolyhedraData(int m, int d)
        {
            int i;
            Polyhedra poly = null;

            poly = new Polyhedra();
            poly.child = null; /* this links the homogenized cone data */
            poly.m = m;
            poly.d = d;
            poly.n = -1;  /* the size of output is not known */
            poly.m_alloc = m + 2; /* the allocated row size of matrix A */
            poly.d_alloc = d;   /* the allocated col size of matrix A */
            poly.numbtype = NumberType.Real;
            InitializeAmatrix(poly.m_alloc, poly.d_alloc, out poly.A);
            InitializeArow(d, out poly.c);           /* cost vector */
            poly.representation = RepresentationType.Inequality;
            poly.homogeneous = false;

            poly.EqualityIndex = new int[m + 2];
            /* size increased to m+2 in 092b because it is used by the child cone, 
               This is a bug fix suggested by Thao Dang. */
            /* ith component is 1 if it is equality, -1 if it is strict inequality, 0 otherwise. */
            for (i = 0; i <= m + 1; i++) poly.EqualityIndex[i] = 0;

            poly.IsEmpty = -1;  /* initially set to -1, neither true nor false, meaning unknown */

            poly.NondegAssumed = false;
            poly.InitBasisAtBottom = false;
            poly.RestrictedEnumeration = false;
            poly.RelaxedEnumeration = false;

            poly.AincGenerated = false;  /* Ainc is a set array to store the input incidence. */

            return poly;
        }

        public bool InitializeConeData(int m, int d, out Cone cone)
        {
            bool success = true;
            int j;

            cone = new Cone();

            /* INPUT: A given representation of a cone: inequality */
            cone.m = m;
            cone.d = d;
            cone.m_alloc = m + 2; /* allocated row size of matrix A */
            cone.d_alloc = d;   /* allocated col size of matrix A, B and Bsave */
            cone.numbtype = NumberType.Real;
            cone.parent = null;

            /* CONTROL: variables to control computation */
            cone.Iteration = 0;

            cone.HalfspaceOrder = RowOrderType.LexMin;

            cone.ArtificialRay = null;
            cone.FirstRay = null;
            cone.LastRay = null; /* The second description: Generator */
            cone.PosHead = null;
            cone.ZeroHead = null;
            cone.NegHead = null;
            cone.PosLast = null;
            cone.ZeroLast = null;
            cone.NegLast = null;
            cone.RecomputeRowOrder = true;
            cone.PreOrderedRun = false;
            cone.GroundSet = new BitArray(cone.m_alloc);
            cone.EqualitySet = new BitArray(cone.m_alloc);
            cone.NonequalitySet = new BitArray(cone.m_alloc);
            cone.AddedHalfspaces = new BitArray(cone.m_alloc);
            cone.WeaklyAddedHalfspaces = new BitArray(cone.m_alloc);
            cone.InitialHalfspaces = new BitArray(cone.m_alloc);
            cone.RayCount = 0;
            cone.FeasibleRayCount = 0;
            cone.WeaklyFeasibleRayCount = 0;
            cone.TotalRayCount = 0;
            cone.ZeroRayCount = 0;
            cone.EdgeCount = 0;
            cone.TotalEdgeCount = 0;
            cone.count_int = 0;
            cone.count_int_good = 0;
            cone.count_int_bad = 0;
            cone.rseed = 1;  /* random seed for random row permutation */

            InitializeBmatrix(cone.d_alloc, out cone.B);
            InitializeBmatrix(cone.d_alloc, out cone.Bsave);
            InitializeAmatrix(cone.m_alloc, cone.d_alloc, out cone.A);

            cone.Edges = new Adjacency[cone.m_alloc];
            cone.InitialRayIndex = new int[d + 1];
            cone.OrderVector = new int[cone.m_alloc + 1];

            cone.newcol = new int[cone.d + 1];
            for (j = 0; j <= cone.d; j++) cone.newcol[j] = j;  /* identity map, initially */
            cone.LinearityDim = -2; /* -2 if it is not computed */
            cone.ColReduced = false;
            cone.d_orig = d;

            /* STATES: variables to represent current state. */
            /*cone.Error;
              cone.CompStatus;
              cone.starttime;
              cone.endtime;
            */

            return success;
        }

        public Cone ConeDataLoad(Polyhedra poly)
        {
            Cone cone = null;
            int d, j;
            int m, i;

            m = poly.m;
            d = poly.d;
            if (!(poly.homogeneous) && poly.representation == RepresentationType.Inequality)
            {
                m = poly.m + 1;
            }
            poly.m1 = m;

            InitializeConeData(m, d, out cone);
            cone.representation = poly.representation;

            /* Points to the original polyhedra data, and reversely */
            cone.parent = poly;
            poly.child = cone;

            for (i = 1; i <= poly.m; i++)
                for (j = 1; j <= cone.d; j++)
                    cone.A[i - 1][j - 1] = poly.A[i - 1][j - 1];

            if (poly.representation == RepresentationType.Inequality && !(poly.homogeneous))
            {
                cone.A[m - 1][0] = 1.0;
                for (j = 2; j <= d; j++) cone.A[m - 1][j - 1] = 0.0;
            }

            return cone;
        }

#if nemkell
void SetLinearity(Matrix M, char *line)
{
  int i=0;
  int eqsize,var;
  char *next;
  const char ct[]=", ";  /* allows separators "," and " ". */

  next=strtok(line,ct);
  eqsize=atol(next); 
  while (i < eqsize && (next=strtok(null,ct))!=null) {
     var=atol(next);
     set_addelem(M.linset,var); i++;
  }
  if (i!=eqsize) {
    fprintf(stderr,"* Warning: there are inconsistencies in linearity setting.\n");
  }
  return;
}

Matrix PolyFile2Matrix (FILE *f, ErrorType *Error)
{
  Matrix M=null;
  int m_input,i;
  int d_input,j;
  RepresentationType rep=Inequality;
  double value;
  bool found=false, newformat=false, successful=false, linearity=false;
  char command[linelenmax], comsave[linelenmax], numbtype[wordlenmax];
  NumberType NT;
//#if !defined(GMPRATIONAL)
  double rvalue;
//#endif

  init(value);
  (*Error)=NoError;
  while (!found)
  {
    if (fscanf(f,"%s",command)==EOF) {
      (*Error)=ImproperInputFormat;
      goto _L99;
    }
    else {
      if (strncmp(command, "V-representation", 16)==0) {
        rep=Generator; newformat=true;
      }
      if (strncmp(command, "H-representation", 16)==0){
        rep=Inequality; newformat=true;
      }
      if (strncmp(command, "partial_enum", 12)==0 || 
          strncmp(command, "equality", 8)==0  ||
          strncmp(command, "linearity", 9)==0 ) {
        linearity=true;
        fgets(comsave,linelenmax,f);
      }
      if (strncmp(command, "begin", 5)==0) found=true;
    }
  }
  fscanf(f, "%ld %ld %s", &m_input, &d_input, numbtype);
  fprintf(stderr,"size = %ld x %ld\nNumber Type = %s\n", m_input, d_input, numbtype);
  NT=GetNumberType(numbtype);
  if (NT==Unknown) {
      (*Error)=ImproperInputFormat;
      goto _L99;
    } 
  M=CreateMatrix(m_input, d_input);
  M.representation=rep;
  M.numbtype=NT;

  for (i = 1; i <= m_input; i++) {
    for (j = 1; j <= d_input; j++) {
      if (NT==Real) {
//#if defined GMPRATIONAL
        *Error=NoRealNumberSupport;
        goto _L99;
//#else
        fscanf(f, "%lf", &rvalue);
        set_d(value, rvalue);
//#endif
      } else {
        fread_rational_value (f, value);
      }
      set(M.matrix[i-1][j - 1],value);
      if (debug) {fprintf(stderr,"a(%3ld,%5ld) = ",i,j); WriteNumber(stderr,value);}
    }  /*of j*/
  }  /*of i*/
  if (fscanf(f,"%s",command)==EOF) {
   	 (*Error)=ImproperInputFormat;
  	 goto _L99;
  }
  else if (strncmp(command, "end", 3)!=0) {
     if (debug) fprintf(stderr,"'end' missing or illegal extra data: %s\n",command);
     (*Error)=ImproperInputFormat;
     goto _L99;
  }
  
  successful=true;
  if (linearity) {
    SetLinearity(M,comsave);
  }
  while (!feof(f)) {
    fscanf(f,"%s", command);
    ProcessCommandLine(f, M, command);
    fgets(command,linelenmax,f); /* skip the CR/LF */
  } 

_L99: ;
  clear(value);
  /* if (f!=null) fclose(f); */
  return M;
}

#endif

        public Polyhedra DDMatrix2Poly(Matrix M, out ErrorType err)
        {
            int i;
            int j;
            Polyhedra poly = null;

            err = ErrorType.NoError;
            if (M.rowsize < 0 || M.colsize < 0)
            {
                err = ErrorType.NegativeMatrixSize;
                return null;
            }
            poly = CreatePolyhedraData(M.rowsize, M.colsize);
            poly.representation = M.representation;
            poly.homogeneous = true;

            for (i = 1; i <= M.rowsize; i++)
            {
                if (M.linset.Get(i))
                {
                    poly.EqualityIndex[i] = 1;
                }
                for (j = 1; j <= M.colsize; j++)
                {
                    poly.A[i - 1][j - 1] = M.matrix[i - 1][j - 1];
                    if (j == 1 && M.matrix[i - 1][j - 1] != 0.0) poly.homogeneous = false;
                }  /*of j*/
            }  /*of i*/
            DoubleDescription(poly, out err);

            return poly;
        }

        public Polyhedra DDMatrix2Poly2(Matrix M, RowOrderType horder, out ErrorType err)
        {
            int i;
            int j;
            Polyhedra poly = null;

            err = ErrorType.NoError;
            if (M.rowsize < 0 || M.colsize < 0)
            {
                err = ErrorType.NegativeMatrixSize;
                return null;
            }
            poly = CreatePolyhedraData(M.rowsize, M.colsize);
            poly.representation = M.representation;
            poly.homogeneous = true;

            for (i = 1; i <= M.rowsize; i++)
            {
                if (M.linset.Get(i))
                {
                    poly.EqualityIndex[i] = 1;
                }
                for (j = 1; j <= M.colsize; j++)
                {
                    poly.A[i - 1][j - 1] = M.matrix[i - 1][j - 1];
                    if (j == 1 && M.matrix[i - 1][j - 1] != 0.0) poly.homogeneous = false;
                }  /*of j*/
            }  /*of i*/
            DoubleDescription2(poly, horder, out err);

            return poly;
        }

        void MatrixIntegerFilter(Matrix M)
        {   /* setting an almost integer to the integer. */
            int i;
            int j;
            double x;

            x = 0.0;
            for (i = 0; i < M.rowsize; i++)
                for (j = 0; j < M.colsize; j++)
                {
                    SnapToInteger(out x, M.matrix[i][j]);
                    M.matrix[i][j] = x;
                }

        }

        void CopyRay(double[] a, int d_origsize, Ray RR,
          RepresentationType rep, int[] reducedcol)
        {
            int j, j1;
            double b;

            b = 0.0;
            for (j = 1; j <= d_origsize; j++)
            {
                j1 = reducedcol[j];
                if (j1 > 0)
                {
                    a[j - 1] = RR.RRay[j1 - 1];
                    /* the original column j is mapped to j1, and thus
                       copy the corresponding component */
                }
                else
                {
                    a[j - 1] = 0.0;
                    /* original column is redundant and removed for computation */
                }
            }

            b = a[0];
            if (rep == RepresentationType.Generator && b != 0)
            {
                a[0] = 1.0;
                for (j = 2; j <= d_origsize; j++)
                    a[j - 1] /= b;    /* normalization for generators */
            }
        }


        void WriteRay(TextWriter f, int d_origsize, Ray RR, RepresentationType rep, int[] reducedcol)
        {
            int j;
            double[] a;


            InitializeArow(d_origsize + 1, out a);

            CopyRay(a, d_origsize, RR, rep, reducedcol);
            for (j = 0; j < d_origsize; j++) WriteNumber(f, a[j]);
            f.WriteLine();
        }

        void WriteArow(TextWriter f, double[] a, int d)
        {
            int j;

            for (j = 0; j < d; j++) WriteNumber(f, a[j]);
            f.WriteLine();
        }

        void WriteAmatrix(TextWriter f, double[][] A, int rowmax, int colmax)
        {
            int i, j;

            if (A == null)
            {
                f.WriteLine("WriteAmatrix: The requested matrix is empty");
                return;
            }
            f.WriteLine("begin");
#if GMPRATIONAL
  fprintf(f, " %ld %ld rational\n",rowmax, colmax);
#else
            f.WriteLine(" {0} {1} real", rowmax, colmax);
#endif
            for (i = 1; i <= rowmax; i++)
            {
                for (j = 1; j <= colmax; j++)
                {
                    WriteNumber(f, A[i - 1][j - 1]);
                }
                f.WriteLine();
            }
            f.WriteLine("end");

        }

        void WriteBmatrix(TextWriter f, int d_size, double[][] B)
        {
            int j1, j2;

            if (B == null)
            {
                f.WriteLine("WriteBmatrix: The requested matrix is empty");
                return;
            }
            for (j1 = 0; j1 < d_size; j1++)
            {
                for (j2 = 0; j2 < d_size; j2++)
                {
                    WriteNumber(f, B[j1][j2]);
                }  /*of j2*/
                f.WriteLine();
            }  /*of j1*/
            f.WriteLine();

        }
#if nemkell
void WriteSetFamily(TextWriter f, SetFamily F)
{
  int i;

  if (F==null){
    f.WriteLine("WriteSetFamily: The requested family is empty");
    return;
  }
  f.WriteLine("begin");
  f.WriteLine("  {0}    {1}", F.famsize, F.setsize);
  for (i=0; i<F.famsize; i++) {
    f.Write(" {0} {1} : ", i+1, F.set[i].Count);
    //*****set_fwrite(f, F.set[i]);
  }
  f.WriteLine("end");

}

void WriteSetFamilyCompressed(FILE *f, SetFamily F)
{
  int i,card;

  if (F==null){
    fprintf(f, "WriteSetFamily: The requested family is empty\n");
    goto _L99;
  }
  fprintf(f,"begin\n");
  fprintf(f,"  %ld    %ld\n", F.famsize, F.setsize);
  for (i=0; i<F.famsize; i++) {
    card=set_card(F.set[i]);
    if (F.setsize - card >= card){
      fprintf(f, " %ld %ld : ", i+1, card);
      set_fwrite(f, F.set[i]);
    } else {
      fprintf(f, " %ld %ld : ", i+1, -card);
      set_fwrite_compl(f, F.set[i]);
    }
  }
  fprintf(f,"end\n");
_L99:;
}
#endif
        void WriteMatrix(TextWriter f, Matrix M)
        {
            int i, linsize;

            if (M == null)
            {
                f.WriteLine("WriteMatrix: The requested matrix is empty");
                return;
            }
            switch (M.representation)
            {
                case RepresentationType.Inequality:
                    f.WriteLine("H-representation"); break;
                case RepresentationType.Generator:
                    f.WriteLine("V-representation"); break;
                case RepresentationType.Unspecified:
                    break;
            }
            linsize = M.linset.Count;
            if (linsize > 0)
            {
                f.Write("linearity {0} ", linsize);
                for (i = 1; i <= M.rowsize; i++)
                    if (M.linset.Get(i)) f.Write(" {0}", i);
                f.WriteLine();
            }
            WriteAmatrix(f, M.matrix, M.rowsize, M.colsize);
            if (M.objective != LPObjectiveType.LPnone)
            {
                if (M.objective == LPObjectiveType.LPmax)
                    f.WriteLine("maximize");
                else
                    f.WriteLine("minimize");
                WriteArow(f, M.rowvec, M.colsize);
            }

        }

        void WriteLP(TextWriter f, LP lp)
        {
            if (lp == null)
            {
                f.WriteLine("WriteLP: The requested lp is empty");
                return;
            }
            f.WriteLine("H-representation");

            WriteAmatrix(f, lp.A, (lp.m) - 1, lp.d);
            if (lp.objective != LPObjectiveType.LPnone)
            {
                if (lp.objective == LPObjectiveType.LPmax)
                    f.WriteLine("maximize");
                else
                    f.WriteLine("minimize");
                WriteArow(f, lp.A[lp.objrow - 1], lp.d);
            }

        }


        void SnapToInteger(out double y, double x)
        {
            /* this is broken.  It does nothing. */
            y = x;
        }


        void WriteReal(TextWriter f, double x)
        {
            /*
          int ix1,ix2,ix;
          double ax;

          ax=x;  
          ix1= (int) (Math.Abs(ax) * 10000.0 + 0.5);
          ix2= (int) (Math.Abs(ax) + 0.5);
          ix2= ix2*10000;
          if ( ix1 == ix2) {
            if (x > 0) {
              ix = (int)(ax + 0.5);
            } else {
              ix = (int)(-ax + 0.5);
              ix = -ix;
            }
            f.Write(" {0}", ix);    //%2ld
          } else
            f.Write(" {0}.9E",ax);   //% .9E
             * */

            f.Write(x);
        }

        void WriteNumber(TextWriter f, double x)
        {
#if GMPRATIONAL
  mpz_t zn,zd;

  mpz_init(zn); mpz_init(zd);
  mpq_canonicalize(x);
  mpq_get_num(zn,x);
  mpq_get_den(zd,x);
  fprintf(f," ");
  if (mpz_sgn(zn)==0){
    fprintf(f,"0");
  } else if (mpz_cmp_ui(zd,1U)==0){
    mpz_out_str(f,10,zn);
  } else {
    mpz_out_str(f,10,zn);fprintf(f,"/");mpz_out_str(f,10,zd);
  }
  mpz_clear(zn); mpz_clear(zd);
#else
            WriteReal(f, x);
#endif
        }

#if nemkell
void WriteIncidence(FILE *f, Polyhedra poly)
{
  SetFamily I;

  switch (poly.representation) {
    case Inequality:
       fprintf(f, "ecd_file: Incidence of generators and inequalities\n");
      break;
    case Generator:
       fprintf(f, "icd_file: Incidence of inequalities and generators\n");
      break;

    default:
      break;
  }
  I=CopyIncidence(poly);
  WriteSetFamilyCompressed(f,I);
  FreeSetFamily(I);
}

void WriteAdjacency(FILE *f, Polyhedra poly)
{
  SetFamily A;

  switch (poly.representation) {
    case Inequality:
       fprintf(f, "ead_file: Adjacency of generators\n");
      break;
    case Generator:
       fprintf(f, "iad_file: Adjacency of inequalities\n");
      break;

    default:
      break;
  }
  A=CopyAdjacency(poly);
  WriteSetFamilyCompressed(f,A);
  FreeSetFamily(A);
}
#endif

        void ComputeAinc(Polyhedra poly)
        {
            /* This generates the input incidence array poly.Ainc, and
               two sets: poly.Ared, poly.Adom. 
            */
            int k;
            int i, m1;
            int j;
            bool redundant;
            Matrix M = null;
            double sum, temp;

            sum = 0.0; temp = 0.0;
            if (poly.AincGenerated == true) return;

            M = CopyOutput(poly);
            poly.n = M.rowsize;
            m1 = poly.m1;
            /* this number is same as poly.m, except when
               poly is given by nonhomogeneous inequalty:
               !(poly.homogeneous) && poly.representation==Inequality,
               it is poly.m+1.   See ConeDataLoad.
            */
            poly.Ainc = new BitArray[m1];
            for (i = 1; i <= m1; i++) poly.Ainc[i - 1] = new BitArray(poly.n);
            poly.Ared = new BitArray(m1);
            poly.Adom = new BitArray(m1);

            for (k = 1; k <= poly.n; k++)
            {
                for (i = 1; i <= poly.m; i++)
                {
                    sum = 0.0;
                    for (j = 1; j <= poly.d; j++)
                    {
                        temp = poly.A[i - 1][j - 1] * M.matrix[k - 1][j - 1];
                        sum += temp;
                    }
                    if (sum == 0.0)
                    {
                        poly.Ainc[i - 1].Set(k, true);
                    }
                }
                if (!(poly.homogeneous) && poly.representation == RepresentationType.Inequality)
                {
                    if (M.matrix[k - 1][0] == 0.0)
                    {
                        poly.Ainc[m1 - 1].Set(k, true);  /* added infinity inequality (1,0,0,...,0) */
                    }
                }
            }

            for (i = 1; i <= m1; i++)
            {
                if (poly.Ainc[i - 1].Count == M.rowsize)
                {
                    poly.Adom.Set(i, true);
                }
            }
            for (i = m1; i >= 1; i--)
            {
                if (poly.Ainc[i - 1].Count == 0)
                {
                    redundant = true;
                    poly.Ared.Set(i, true);
                }
                else
                {
                    redundant = false;
                    for (k = 1; k <= m1; k++)
                    {
                        if (k != i && !poly.Ared.Get(k) && !poly.Adom.Get(k) &&
                            poly.Ainc[i - 1].IsSubsetOf(poly.Ainc[k - 1]))
                        {    //****
                            if (!redundant)
                            {
                                redundant = true;
                            }
                            poly.Ared.Set(i, true);
                        }
                    }
                }
            }

            poly.AincGenerated = true;


        }

        bool InputAdjacentQ(Polyhedra poly,
          int i1, int i2)
        /* Before calling this function, RedundantSet must be 
           a set of row indices whose removal results in a minimal
           nonredundant system to represent the input polyhedron,
           DominantSet must be the set of row indices which are
           active at every extreme points/rays.
        */
        {
            bool adj = true;
            int i;
            BitArray common;
            int lastn = 0;

            if (poly.AincGenerated == false) ComputeAinc(poly);
            if (lastn != poly.n)
            {
                //if (lastn >0) set_free(common);
                common = new BitArray(poly.n);
                lastn = poly.n;
            }
            if (poly.Ared.Get(i1) || poly.Ared.Get(i2))
            {
                adj = false;
                return adj;
            }
            if (poly.Adom.Get(i1) || poly.Adom.Get(i2))
            {
                // dominant inequality is considered adjacencent to all others.
                adj = true;
                return adj;
            }
            common = poly.Ainc[i1 - 1].And(poly.Ainc[i2 - 1]);       //**** intersection
            i = 0;
            while (i < poly.m1 && adj == true)
            {
                i++;
                if (i != i1 && i != i2 && !poly.Ared.Get(i) &&
                    !poly.Adom.Get(i) && common.IsSubsetOf(poly.Ainc[i - 1]))
                {
                    adj = false;
                }
            }

            return adj;
        }

#if nemkell
void WriteInputIncidence(FILE *f, Polyhedra poly)
{
  SetFamily I;

  if (poly.AincGenerated==false) ComputeAinc(poly);
  switch (poly.representation) {
  case Inequality:
    fprintf(f,"icd_file: Incidence of inequalities and generators\n");
    break;

  case Generator:
   fprintf(f,"ecd_file: Incidence of generators and inequalities\n");
    break;

  default:
    break;
  }

  I=CopyInputIncidence(poly);
  WriteSetFamilyCompressed(f,I);
  FreeSetFamily(I);
}

void WriteInputAdjacency(FILE *f, Polyhedra poly)
{
  SetFamily A;

  if (poly.AincGenerated==false){
    ComputeAinc(poly);
  }
  switch (poly.representation) {
  case Inequality:
    fprintf(f, "iad_file: Adjacency of inequalities\n");
    break;

  case Generator:
    fprintf(f, "ead_file: Adjacency of generators\n");
    break;

  default:
    break;
  }
  A=CopyInputAdjacency(poly);
  WriteSetFamilyCompressed(f,A);
  FreeSetFamily(A);
}


void WriteProgramDescription(FILE *f)
{
  fprintf(f, "* cddlib: a double description library:%s\n", DDVERSION);
  fprintf(f, "* compiled for %s arithmetic.\n", ARITHMETIC);
  fprintf(f,"* %s\n",COPYRIGHT);
}

void WriteTimes(FILE *f, time_t starttime, time_t endtime)
{
  int ptime,ptime_sec,ptime_minu, ptime_hour;
  
/* 
   ptime=difftime(endtime,starttime); 
   This function is ANSI standard, but not available sometime 
*/
  ptime=(endtime - starttime);     
 /* This is to replace the line above, but it may not give 
    correct time in seconds */ 
  ptime_hour=ptime/3600;
  ptime_minu=(ptime-ptime_hour*3600)/60;
  ptime_sec=ptime%60;
  fprintf(f,"* Computation started at %s",asctime(localtime(&starttime)));
  fprintf(f,"*             ended   at %s",asctime(localtime(&endtime)));
  fprintf(f,"* Total processor time = %ld seconds\n",ptime);
  fprintf(f,"*                      = %ld h %ld m %ld s\n",ptime_hour, ptime_minu, ptime_sec);
}

void WriteDDTimes(FILE *f, Polyhedra poly)
{
  WriteTimes(f,poly.child.starttime,poly.child.endtime);
}

void WriteLPTimes(FILE *f, LPPtr lp)
{
  WriteTimes(f,lp.starttime,lp.endtime);
}

void WriteLPStats(FILE *f)
{
  time_t currenttime;
  
  time(&currenttime);
  
  fprintf(f, "\n*--- Statistics of pivots ---\n");
//#if defined GMPRATIONAL
  fprintf(f, "* f0 = %ld (float basis finding pivots)\n",ddf_statBApivots);
  fprintf(f, "* fc = %ld (float CC pivots)\n",ddf_statCCpivots);
  fprintf(f, "* f1 = %ld (float dual simplex phase I pivots)\n",ddf_statDS1pivots);
  fprintf(f, "* f2 = %ld (float dual simplex phase II pivots)\n",ddf_statDS2pivots);
  fprintf(f, "* f3 = %ld (float anticycling CC pivots)\n",ddf_statACpivots);
  fprintf(f, "* e0 = %ld (exact basis finding pivots)\n",statBApivots);
  fprintf(f, "* ec = %ld (exact CC pivots)\n",statCCpivots);
  fprintf(f, "* e1 = %ld (exact dual simplex phase I pivots)\n",statDS1pivots);
  fprintf(f, "* e2 = %ld (exact dual simplex phase II pivots)\n",statDS2pivots);
  fprintf(f, "* e3 = %ld (exact anticycling CC pivots)\n",statACpivots);
  fprintf(f, "* e4 = %ld (exact basis verification pivots)\n",statBSpivots);
//#else
  fprintf(f, "f0 = %ld (float basis finding pivots)\n",statBApivots);
  fprintf(f, "fc = %ld (float CC pivots)\n",statCCpivots);
  fprintf(f, "f1 = %ld (float dual simplex phase I pivots)\n",statDS1pivots);
  fprintf(f, "f2 = %ld (float dual simplex phase II pivots)\n",statDS2pivots);
  fprintf(f, "f3 = %ld (float anticycling CC pivots)\n",statACpivots);
//#endif
 WriteLPMode(f);
  WriteTimes(f,statStartTime, currenttime);
}

void WriteLPMode(FILE *f)
{
  fprintf(f, "\n* LP solver: ");
  switch (choiceLPSolverDefault) {
    case DualSimplex:
      fprintf(f, "DualSimplex\n");
      break;
    case CrissCross:
      fprintf(f, "Criss-Cross\n");
      break;
    default: break;
  }
  
  fprintf(f, "* Redundancy cheking solver: ");
  switch (choiceRedcheckAlgorithm) {
    case DualSimplex:
      fprintf(f, "DualSimplex\n");
      break;
    case CrissCross:
      fprintf(f, "Criss-Cross\n");
      break;
    default: break;
  }
  
  fprintf(f, "* Lexicographic pivot: ");
  if (choiceLexicoPivotQ)  fprintf(f, " on\n"); 
  else fprintf(f, " off\n"); 

}


void WriteRunningMode(FILE *f, Polyhedra poly)
{
  if (poly.child!=null){
    fprintf(f,"* roworder: ");
    switch (poly.child.HalfspaceOrder) {

    case MinIndex:
      fprintf(f, "minindex\n");
      break;

    case MaxIndex:
      fprintf(f, "maxindex\n");
      break;

    case MinCutoff:
      fprintf(f, "mincutoff\n");
      break;

    case MaxCutoff:
      fprintf(f, "maxcutoff\n");
    break;

    case MixCutoff:
      fprintf(f, "mixcutoff\n");
      break;

    case LexMin:
      fprintf(f, "lexmin\n");
      break;

    case LexMax:
      fprintf(f, "lexmax\n");
      break;

    case RandomRow:
      fprintf(f, "random  %d\n",poly.child.rseed);
      break;

    default: break;
    }
  }
}


void WriteCompletionStatus(FILE *f, Cone cone)
{
  if (cone.Iteration<cone.m && cone.CompStatus==AllFound) {
    fprintf(f,"*Computation completed at Iteration %4ld.\n", cone.Iteration);
  } 
  if (cone.CompStatus == RegionEmpty) {
    fprintf(f,"*Computation completed at Iteration %4ld because the region found empty.\n",cone.Iteration);
  }   
}

void WritePolyFile(FILE *f, Polyhedra poly)
{
  WriteAmatrix(f,poly.A,poly.m,poly.d);
}


void WriteErrorMessages(FILE *f, ErrorType Error)
{
  switch (Error) {

  case DimensionTooLarge:
    fprintf(f, "*Input Error: Input matrix is too large:\n");
    fprintf(f, "*Please increase MMAX and/or NMAX in the source code and recompile.\n");
    break;

  case IFileNotFound:
    fprintf(f, "*Input Error: Specified input file does not exist.\n");
    break;

  case OFileNotOpen:
    fprintf(f, "*Output Error: Specified output file cannot be opened.\n");
    break;

  case NegativeMatrixSize:
    fprintf(f, "*Input Error: Input matrix has a negative size:\n");
    fprintf(f, "*Please check rowsize or colsize.\n");
    break;

  case ImproperInputFormat:
    fprintf(f,"*Input Error: Input format is not correct.\n");
    fprintf(f,"*Format:\n");
    fprintf(f," begin\n");
    fprintf(f,"   m   n  NumberType(real, rational or integer)\n");
    fprintf(f,"   b  -A\n");
    fprintf(f," end\n");
    break;

  case EmptyVrepresentation:
    fprintf(f, "*Input Error: V-representation is empty:\n");
    fprintf(f, "*cddlib does not accept this trivial case for which output can be any inconsistent system.\n");
    break;

  case EmptyHrepresentation:
    fprintf(f, "*Input Error: H-representation is empty.\n");
    break;

  case EmptyRepresentation:
    fprintf(f, "*Input Error: Representation is empty.\n");
    break;

  case NoLPObjective:
    fprintf(f, "*LP Error: No LP objective (max or min) is set.\n");
    break;

  case NoRealNumberSupport:
    fprintf(f, "*LP Error: The binary (with GMP Rational) does not support Real number input.\n");
    fprintf(f, "         : Use a binary compiled without -DGMPRATIONAL option.\n");
    break;

 case NotAvailForH:
    fprintf(f, "*Error: A function is called with H-rep which does not support an H-representation.\n");
    break;

 case NotAvailForV:
    fprintf(f, "*Error: A function is called with V-rep which does not support an V-representation.\n");
    break;

 case CannotHandleLinearity:
    fprintf(f, "*Error: The function called cannot handle linearity.\n");
    break;

 case RowIndexOutOfRange:
    fprintf(f, "*Error: Specified row index is out of range\n");
    break;

 case ColIndexOutOfRange:
    fprintf(f, "*Error: Specified column index is out of range\n");
    break;

 case LPCycling:
    fprintf(f, "*Error: Possibly an LP cycling occurs.  Use the Criss-Cross method.\n");
    break;
    
 case NumericallyInconsistent:
    fprintf(f, "*Error: Numerical inconsistency is found.  Use the GMP exact arithmetic.\n");
    break;
    
  case NoError:
    fprintf(f,"*No Error found.\n");
    break;
  }
}
#endif

        SetFamily CopyIncidence(Polyhedra poly)
        {
            SetFamily F = null;
            int k;
            int i;

            if (poly.child == null || poly.child.CompStatus != CompStatusType.AllFound) return F;
            if (poly.AincGenerated == false) ComputeAinc(poly);
            F = CreateSetFamily(poly.n, poly.m1);
            for (i = 1; i <= poly.m1; i++)
                for (k = 1; k <= poly.n; k++)
                    if (poly.Ainc[i - 1].Get(k)) F.set[k - 1].Set(i, true);

            return F;
        }

        SetFamily CopyInputIncidence(Polyhedra poly)
        {
            int i;
            SetFamily F = null;

            if (poly.child == null || poly.child.CompStatus != CompStatusType.AllFound) return F;
            if (poly.AincGenerated == false) ComputeAinc(poly);
            F = CreateSetFamily(poly.m1, poly.n);
            for (i = 0; i < poly.m1; i++)
            {
                poly.Ainc[i].CopyTo(F.set[i]);
            }

            return F;
        }

        SetFamily CopyAdjacency(Polyhedra poly)
        {
            Ray RayPtr1, RayPtr2;
            SetFamily F = null;
            int pos1, pos2;
            int lstart, k, n;
            BitArray linset, allset;
            bool adj;

            if (poly.n == 0 && poly.homogeneous && poly.representation == RepresentationType.Inequality)
            {
                n = 1; /* the origin (the unique vertex) should be output. */
            }
            else n = poly.n;
            linset = new BitArray(n);
            allset = new BitArray(n);
            if (poly.child == null || poly.child.CompStatus != CompStatusType.AllFound) return F;
            F = CreateSetFamily(n, n);
            if (n <= 0) return F;
            poly.child.LastRay.Next = null;
            for (RayPtr1 = poly.child.FirstRay, pos1 = 1; RayPtr1 != null;
                          RayPtr1 = RayPtr1.Next, pos1++)
            {
                for (RayPtr2 = poly.child.FirstRay, pos2 = 1; RayPtr2 != null;
                                RayPtr2 = RayPtr2.Next, pos2++)
                {
                    if (RayPtr1 != RayPtr2)
                    {
                        CheckAdjacency(poly.child, RayPtr1, RayPtr2, out adj);
                        if (adj)
                        {
                            F.set[pos1 - 1].Set(pos2, true);
                        }
                    }
                }
            }
            lstart = poly.n - poly.ldim + 1;
            allset = allset.Not();  /* allset is set to the ground set. */
            for (k = lstart; k <= poly.n; k++)
            {
                linset.Set(k, true);     /* linearity set */
                allset.CopyTo(F.set[k - 1]);  /* linearity generator is adjacent to all */
            }
            for (k = 1; k < lstart; k++)
            {
                F.set[k - 1] = F.set[k - 1].Or(linset); // union
                /* every generator is adjacent to all linearity generators */
            }

            return F;
        }

        SetFamily CopyInputAdjacency(Polyhedra poly)
        {
            int i, j;
            SetFamily F = null;

            if (poly.child == null || poly.child.CompStatus != CompStatusType.AllFound) return F;
            if (poly.AincGenerated == false) ComputeAinc(poly);
            F = CreateSetFamily(poly.m1, poly.m1);
            for (i = 1; i <= poly.m1; i++)
            {
                for (j = 1; j <= poly.m1; j++)
                {
                    if (i != j && InputAdjacentQ(poly, i, j))
                    {
                        F.set[i - 1].Set(j, true);
                    }
                }
            }

            return F;
        }

        Matrix CopyOutput(Polyhedra poly)
        {
            Ray Ray;
            Matrix M = null;
            int i = 0, total;
            int j, j1;
            double b;
            RepresentationType outputrep = RepresentationType.Inequality;
            bool outputorigin = false;

            b = 0.0;
            total = poly.child.LinearityDim + poly.child.FeasibleRayCount;

            if (poly.child.d <= 0 || poly.child.newcol[1] == 0) total = total - 1;
            if (poly.representation == RepresentationType.Inequality) outputrep = RepresentationType.Generator;
            if (total == 0 && poly.homogeneous && poly.representation == RepresentationType.Inequality)
            {
                total = 1;
                outputorigin = true;
                /* the origin (the unique vertex) should be output. */
            }
            if (poly.child == null || poly.child.CompStatus != CompStatusType.AllFound) return M;

            M = CreateMatrix(total, poly.d);
            Ray = poly.child.FirstRay;
            while (Ray != null)
            {
                if (Ray.feasible)
                {
                    CopyRay(M.matrix[i], poly.d, Ray, outputrep, poly.child.newcol);
                    i++;  /* 086 */
                }
                Ray = Ray.Next;
            }
            for (j = 2; j <= poly.d; j++)
            {
                if (poly.child.newcol[j] == 0)
                {
                    /* original column j is dependent on others and removed for the cone */
                    b = poly.child.Bsave[0][j - 1];
                    if (outputrep == RepresentationType.Generator && b > 0.0)
                    {
                        M.matrix[i][0] = 1.0;  /* Normalize */
                        for (j1 = 1; j1 < poly.d; j1++)
                            M.matrix[i][j1] = poly.child.Bsave[j1][j - 1] / b;
                    }
                    else
                    {
                        for (j1 = 0; j1 < poly.d; j1++)
                            M.matrix[i][j1] = poly.child.Bsave[j1][j - 1];
                    }
                    M.linset.Set(i + 1, true);
                    i++;
                }
            }
            if (outputorigin)
            {
                /* output the origin for homogeneous H-polyhedron with no rays. */
                M.matrix[0][0] = 1.0;
                for (j = 1; j < poly.d; j++)
                {
                    M.matrix[0][j] = 0.0;
                }
            }
            MatrixIntegerFilter(M);
            if (poly.representation == RepresentationType.Inequality)
                M.representation = RepresentationType.Generator;
            else
                M.representation = RepresentationType.Inequality;


            return M;
        }

        Matrix CopyInput(Polyhedra poly)
        {
            Matrix M = null;
            int i;

            M = CreateMatrix(poly.m, poly.d);
            CopyAmatrix(M.matrix, poly.A, poly.m, poly.d);

            for (i = 1; i <= poly.m; i++)
                if (poly.EqualityIndex[i] == 1) M.linset.Set(i, true);
            MatrixIntegerFilter(M);
            if (poly.representation == RepresentationType.Generator)
                M.representation = RepresentationType.Generator;
            else
                M.representation = RepresentationType.Inequality;
            return M;
        }

        public Matrix CopyGenerators(Polyhedra poly)
        {
            Matrix M = null;

            if (poly.representation == RepresentationType.Generator)
            {
                M = CopyInput(poly);
            }
            else
            {
                M = CopyOutput(poly);
            }
            return M;
        }

        Matrix CopyInequalities(Polyhedra poly)
        {
            Matrix M = null;

            if (poly.representation == RepresentationType.Inequality)
            {
                M = CopyInput(poly);
            }
            else
            {
                M = CopyOutput(poly);
            }
            return M;
        }

#if nemkell
/****************************************************************************************/
/*  rational number (a/b) read is taken from Vinci by Benno Bueeler and Andreas Enge    */
/****************************************************************************************/
void sread_rational_value (char *s, double value)
   /* reads a rational value from the specified string "s" and assigns it to "value"    */
   
{
   char     *numerator_s=null, *denominator_s=null, *position;
   int      sign = 1;
   double   numerator, denominator;
//#if defined GMPRATIONAL
   mpz_t znum, zden;
//#else
   double rvalue;
//#endif
  
   /* determine the sign of the number */
   numerator_s = s;
   if (s [0] == '-')
   {  sign = -1;
      numerator_s++;
   }
   else if (s [0] == '+')
      numerator_s++;
      
   /* look for a sign '/' and eventually split the number in numerator and denominator */
   position = strchr (numerator_s, '/');
   if (position != null)
   {  *position = '\0'; /* terminates the numerator */
      denominator_s = position + 1;
   };

   /* determine the floating point values of numerator and denominator */
   numerator=atol (numerator_s);
  
   if (position != null)
   { 
     denominator=atol (denominator_s);  
   }
   else denominator = 1;

/* 
   fprintf(stderr,"\nrational_read: numerator %f\n",numerator);
   fprintf(stderr,"rational_read: denominator %f\n",denominator);
   fprintf(stderr,"rational_read: sign %d\n",sign); 
*/

//#if defined GMPRATIONAL
   mpz_init_set_str(znum,numerator_s,10);
   if (sign<0) mpz_neg(znum,znum);
   mpz_init(zden); mpz_set_ui(zden,1);
   if (denominator_s!=null) mpz_init_set_str(zden,denominator_s,10);
   mpq_set_num(value,znum); mpq_set_den(value,zden);
   mpq_canonicalize(value);
   mpz_clear(znum); mpz_clear(zden);
   /* num=(int)sign * numerator; */
   /* den=(unsigned int) denominator; */
   /* mpq_set_si(value, num, den); */
//#elif defined GMPFLOAT
   rvalue=sign * numerator/ (signed int) denominator;
   mpf_set_d(value, rvalue);
//#else
   rvalue=sign * numerator/ (signed int) denominator;
   dset_d(value, rvalue);
//#endif
   if (debug) {
     fprintf(stderr,"rational_read: "); 
     WriteNumber(stderr,value); fprintf(stderr,"\n");
   }
}
   

void fread_rational_value (FILE *f, double value)
   /* reads a rational value from the specified file "f" and assigns it to "value"      */
   
{
   char     number_s [255];
   double rational_value;
   
   init(rational_value);
   fscanf(f, "%s ", number_s);
   sread_rational_value (number_s, rational_value);
   set(value,rational_value);
   clear(rational_value);
}
   
/****************************************************************************************/


/* end of cddio.c */
#endif

    }
}
